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Abstract 

A numerical technique used to solve boundary value problems is modified to find periodic 
steady-state solutions of nonautonomous dynamical systems. The technique uses a matrix 
representation of the time derivative obtained through trigonometric interpolation of a 
periodic function. Such a differentiation matrix yields exact values for the derivative of a 
trigonometric polynomial and therefore, can be used as the main ingredient of a numerical 
method to solve nonlinear dynamical systems. We apply this technique to obtain some 
limit cycles and bifurcation points of a sinusoidally driven pendulum and the steady-state 
response of an electric circuit. 
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1 Introduction 



Dynamical systems can be given by a set of ordinary differential equations and some initial 
conditions. Such problems are usually solved by conventional procedures such as Runge- 
Kutta methods, multistep methods or general linear methods (see for example jT]). For 
instance, shooting methods j2], extrapolation methods |2J and Newton methods based on 
the Poincare map [3] has been applied to accelerate the convergence of the solution to the 
periodic steady-state of electric systems. Since some of these implementations requires 
a considerable computational effort, an alternative and more simple technique may be 
desired. 

In a sequel of papers (see |S] and references therein), a Galerkin-collocation-type method 
for solving differential boundary value problems has been developed. In some restricted 
cases this technique yields exact results, i.e., the numerical output can be interpolated to 
obtain the functions that solve exactly the problem, and therefore it have been used to 
implement a numerical scheme to solve dynamical systems that consists basically in the 
substitution of the derivatives appearing in the differential equations by finite-dimensional 
matrices whose entries depend on certain set of points in a simple form. In some cases, the 
nodes have to be chosen according to some criterion in order to accelerate the convergence 
to the solution. The differentiation matrices used in this technique arise naturally in the 
context of the interpolation of functions. 

The aim of this paper is to present a novel method to obtain numerical approximants to 
the steady-state solutions of nonautonomous systems, based on the use of the differenti- 
ation matrix for periodic functions given in 



The main idea of the method is established in Sec. Eland two examples are given in Sec. H3 



2 Discrete formalism 

In this section we only present the main results of our discretization scheme for periodic 
functions; proofs and further applications can be found in [H] . Let D be the N x N matrix 
whose components are given by 



D jk 




J = k, 




where — ir < ti < t2 < ■ ■ ■ < < it and 
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Let x(t) be a trigonometric polynomial of degree at most n and x and x' denote the 
iV x 1 vectors whose elements are x(tj) and x'(tj), i.e., the values of the polynomial and 
its derivative at the nodes respectively. Then, the matrix D applied to x becomes equal 
to x' whenever N = 2n + 1. Since the derivative of a trigonometric polynomial is again a 
trigonometric polynomial of the same degree, we have that 



x 



{k) = D k x, k = 0,1,2,..., (1) 



where x^ is the vector whose elements are given by the kth. derivative of x(t) evaluated 
at the nodes. 

The differentiation matrix D takes the simple form 

r o, j = k, 

D jk = I (-iy +k . ( 2 ) 

{2sm%(j-ky 
if the nodes are chosen to be the N equidistant points 1 

^ = -7T + j = l,2,---,iV. (3) 

If the periodic function x(t) is not a polynomial, a residual vector depending on x(t), k 
and N must be added to the right-hand side of (JTJ). However, it is expected that the norm 
of such a vector approaches zero as the number of nodes is increased since x(t) can be 
expanded in a Fourier series. Therefore, if N is great enough, the residual vector can be 
ignored and the function x^ k ' (t) can be approximated by an interpolation of the elements 
of D k x. 

Let us consider now a nonautonomous system of m components described by 

x = f(x,wt), (4) 

where x is the vector of components xi(t),X2(t), . . . ,x m (t) and / is a nonlinear vector 
function of m + 1 variables, periodic in t, with period T = 2tc/uj. A very important 
problem is to find the response of the system in the steady-state regime. Taking into 
account the periodicity of this response, we can use the differentiation matrix (J2J) to obtain 
approximants to the steady-state solution of (@J) according to the following scheme. 
First of all, we need to change ut — > t in (jlj in order to change the period of / to 2n. 
Thus, (dJ) becomes 

ux = f(x,t), (5) 

Let us take an odd number of points tj as given by and evaluate (0) at each node 
to form an equality between Nm x 1 vectors in such a way that the first A^ entries of 
the left-hand side are the components of the vector xi, i.e., Xi(ti), Xi(t?), . . . , ii(tjv); the 

1 The same result is obtained for the set of points tj = n(2j — N — 1)/N. 
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following iV entries are the components of x 2 , i.e., x 2 (t 1 ) , x 2 (t 2 ) , • • • ,^2(^)1 and so on. 
Now we can approximate the vectors blocks Xk by Dxk to obtain the discrete form of (J5J), 
which can be written as 

N 

u^DjiXkih) = f k (x 1 (t j ),x 2 , (tj) . . .,x m (tj),tj), j = 1,2, . . .,N, (6) 
1=1 

where k — 1, 2, . . . , to, or in the more compact form 

wDX = F, (7) 

where F denotes the Nm x 1 vector whose elements are given by the right-hand side of 
(JBJ) first running j and then k, X is the vector whose elements are given by Xk(tj) (the 
unknown solution) ordered in a similar way and D = l m Cg) D, where l m is the identity 
matrix of dimension m. 

The points on which our method is based are the following: 

1. If (@J) has a limit cycle, then the solution of (J7J) is an approximation of the steady- 
state solution of (HJ). 

2. Since (J7J) [or equivalently (jUJ)] is a system of Nm nonlinear equations with Nm un- 
knowns Xk(tj), its solution can be obtained by using a standard procedure (Newton's 
method for instance). 

This shows that a nonautonomous system in the steady-state regime can be described 
approximately by a system of nonlinear algebraic equations. The procedure sketched in 
these statements is not concerned at all with the initial conditions of the system and 
yields simultaneously all the values of Xk(tj) (at all times). Thus, this method is quite 
different in essence to those designed as initial-value-problem solvers. 
To obtain the solution of (jUJ) we can use the Newton method or some variation of it. How- 
ever, as it is well-known, not always is easy to give a good initial approximation to attain 
convergence. A good issue is to make a few Runge-Kutta integrations of the system to 
yield an initial approximation Xo. Since this procedure is time-consuming, a more simple 
way to find X should be tried, but it will depend on the problem to solve. 
Due to the nature of our method to approximate the steady-state solution of (j5J), an al- 
gorithm for this technique will consist necessarily in the algorithm selected to solve the 
set of nonlinear equations. 

In the following section we put this method at work on two important nonautonomous 
dynamical systems. 



3 Test cases 

We have chosen a chaotic mechanical system and a electric circuit of actual use as test 
cases. To find their steady-state solutions, we rewrite the equations describing the dy- 
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namics of these problems in the discretized form (JJJ) and use a FORTRAN90 program 
and standard libraries running on a personal computer to solve them. 



3.1 Test case 1: the driven pendulum 

In spite that the classical pendulum is a very old problem, interest on it is still growing. 
The pendulum becomes a chaotic system when it is driven at the pivot point. Let us 
consider a rigid and planar pendulum consisting in a mass attached to a light rod of 
length I which is vertically driven by a sinusoidal force of the form — A cos u>t and damped 
by a linear viscous force with damping p, as in [7j. If 9 denotes the angular displacement 
of the pendulum measured from the vertically downward position, the equation of motion 
is 

■jp +a— + (1 + 6 cos cut) sin = 0, (8) 

where 

2p Au 2 



>/lg' l ' 

To compare the steady-state solutions yielded by our procedure with those obtained by 
other authors we take < b < 200 (the driving amplitude) as the control parameter, 
a = 0.1 and uj = 17.5. For values of b near zero, a good initial guess is a vector X whose 
entries are close to tt. To catch more solutions for a given value of b, one can try a vector 
with a sinusoidal deviation from tt. Once we have found a solution X(p), we proceed to 
obtain X(p + Sp) by using X = X(p). We have taken N = 101 nodes of the form (jSJ) in 
this case. 

As it is known, this system presents period-doubling and bifurcation. Figure 1 displays 
these phenomena for solutions of (JBJ) close to the inverted state 9 = tt. To make the 
drawing more simple, we only plot the maximum and minimum values of the angular 
displacement 9{t) as functions of the control parameter b. Other solutions, as the hanging 
state, are not considered here. 

The displayed information is in line with the findings of other authors (see for example 
but in our case the problem of stability of the solution is not revised. In Fig. 1, 
9\ is the solution corresponding to the inverted state, 9<i and 6*3 correspond in essence to 
the same period- 1 solution (one can get the second from the first through a reflection on 
the axis 9 = tt) and #4 and #5 are period-2 solutions oscillating both of them about the 
vertical. The asymmetrical way in which # 5 oscillates is illustrated in Fig. 2. 



5 



1.8 
1.6 
1.4 
1.2 

Bin 1 

0.8 
0.6 
0.4 
0.2 



I I I I I I I 
mux 

w=17.5 — 

a=0.1 — — " ~~ 

— min 


1 1 

jnax 
V5 




/,max 
/ 04 






min ' — • 

^ — ^^2* 

1 1 1 1 1 1 1 1 



20 40 60 80 100 120 140 160 180 200 



Figure 1: Maximum and minimum values of the angular displacement 9 of the 
vertically driven pendulum © vs. the driven amplitude b for a = 0.1 and uj = 17.5. 
For any value of b, the solution 6\ = it (the inverted pendulum) always was found. 
The curves labeled by 6^, k = 2, • • • , 5 correspond to different solutions. 
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3.2 Test case 2: a commutation circuit 

The actual interest in circuits as the one depicted in Fig. 3, relies on the extended use of 
them in modern power supplies. Usually a primary DC power supply (a battery or the 
rectified and filtered 60 Hz main line) is commutated at high frequencies (tens of KHz 
to MHz) in order to obtain a pulse waveform, which can be coupled to the load in some 
way to obtain a predetermined load voltage. On the circuit shown in Fig. 3, V s is a 
simplified model of a transformer secondary, which is followed by a rectifier-filter circuit 
and a resistive load (R4). R\ stands for the series resistance of the diode and R2 and R3 
model the equivalent series resistance of C\ and C2 respectively. R2 and R3 are important 
parameters because they increase ripple voltages and circuit losses. 
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Figure 3: A typical commutation circuit. 



The equations for the state variables V\ 
the form 



x 2 



£3 



C\Ri 



V s - Xi 



Vd-isR^ 



C2{R?, + Ra) 



R 1 



exp( 

x 2 + RiX 3 j , 

R3R4 



xi, V 2 = x 2 and ii 
V d q, 



x 3 can be simplified to 



r]kT' 



X2 - 



(i?3 + R4) (-R3 + R4 



-X3 



-Vd- i s -Ri(expi 



rjkT' 



- 1 



where 



V d = V s - C x (Ri + R 2 )x\ - xi - Rix 3 , 



(9) 



k is the Boltzmann constant, q is the electron charge, T the absolute temperature which 
is taken at the standard value T = 300°i^ and rj is the emission coefficient. The driving 
voltage is chosen as a square wave pulse, i.e., V s (t) = A m sgn(27rt/T), t G [— T/2,T/2], 
and the remaining parameters of this problem are taken at typical values: A m = 5.6V, 



7 



T = 10" 5 s, i s = 10~ 8 A, R x = 0.0149ft, R 2 = 0.15ft, R 3 = 0.2ft, R 4 = 2.0ft, d = 
470.0 x 10" 6 F, C 2 = 20.0 x 10~ 6 F, L = 20.0 x 10~ 6 H, r] = 0.8953. 

In order to have results to compare with, we simulated the circuit using a LT version of 
the Spice program which is currently available for free from Linear Technology, a major 
semiconductor devices manufacturer. Spice 8 J is a standard tool used for electric circuits 
simulation. Such a program compute the steady-state response in a sequential way and 
therefore, through a number of periods of the excitation. For this circuit, the response of 
the system was computed by the Spice simulation over 150 cycles to achieve the steady- 
state. A maximum step size of 0.04yus was specified. In our case, we use 251 equispaced 
time values tj of the form (J3J) (yielding a constant step size approximately equal to 0.04/xs) 
and all of the steady-state response values were computed at once. Since Spice uses an 
adaptive strategy for the step size, it is not possible to compare both of the results in a 
point-by-point basis. However, it comes out that they agree globally up to 10~ 2 . 
The output yielded by the present method and the Spice simulation for the more used 
variables id = £3 + C\X\ and Vo = C2R3X2 + %2, is displayed in Figs. (4a) and (4b), 
respectively. 




Figure 4: Steady-state solution for the diode current id = £3 + C\x\ and the voltage 
Vo = C2R3X2 + X2 along two periods obtained by the present method with N = 251 
nodes (4a) and by a Spice simulation (4b). The left-hand side and right-hand side 
vertical axes are used for different scales. 



4 Final remarks 



As the simple tests of the previous section show, the method presented in this paper gives 
a novel approach to obtain the limit cycles of nonautonomous systems and represents an 
alternative to the conventional methods of integration. Based on a well-behaved differen- 
tiation matrix for periodic functions, it departs from the standard procedures in the way 
in which the limit cycle is found: the usual integration of the system for long times is 



8 



replaced by the problem of finding the solution of a set of nonlinear algebraic equations. 
This feature makes this method a tool for studying dynamical systems from a new point 
of view. 
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